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ABSTRACT 

Cosmological simulations make use of sub-grid recipes for the implementation of galac- 
tic winds driven by massive stars because direct injection of supernova energy in ther- 
mal form leads to strong radiative losses, rendering the feedback inefficient. We argue 
that the main cause of the catastrophic cooling is a mismatch between the mass of 
the gas in which the energy is injected and the mass of the parent stellar population. 
Because too much mass is heated, the temperatures are too low and the cooling times 
too short. We use analytic arguments to estimate, as a function of the gas density 
and the numerical resolution, the minimum heating temperature that is required for 
the injected thermal energy to be efficiently converted into kinetic energy. We then 
propose and test a stochastic implementation of thermal feedback that uses this min- 
imum temperature increase as an input parameter and that can be employed in both 
particle- and grid-based codes. We use smoothed particle hydrodynamics simulations 
to test the method on models of isolated disc galaxies in dark matter haloes with 
total mass 10^^ and 10^^ h~^ Mq. The thermal feedback strongly suppresses the star 
formation rate and can drive massive, large-scale outflows without the need to turn off 
radiative cooling temporarily. In accord with expectations derived from analytic ar- 
guments, for sufficiently high resolution the results become insensitive to the imposed 
temperature jump and also agree with high-resolution simulations employing kinetic 
feedback. 

Key words: methods: numerical — ISM: bubbles — ISM: jets and outflows — galax- 
ies: evolution — galaxies: formation — galaxies: ISM 
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1 INTRODUCTION 

It is widely accepted that star formation (SF) feeds back 
energy into the interstellar medium (ISM). The energy re- 
leased by massive stars, both through stellar winds and core- 
collapse supernova (SN) explosions, can efficiently suppress 
SF by evaporating dense, star-forming clouds, by generating 
supersonic turbulence and, eventually, by generating pow- 
erful, large-scale outflows that eject gas from galaxies and 
enrich the intergalactic medium. 

Modern cosmological simulations that follow the for- 
mation and evolution of galaxies still lack both the resolu- 
tion and the physics that is required to model the multi- 
phase ISM and individual massive stars or SN explosions. 
The same is true for simulations of individual galaxies, al- 
though the resolution of such models is now sufficient to 
begin to crudely disentangle the relative roles of the differ- 
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ent mechanisms through which massive stars inject energy 
and momentum (e.g. radiation pres sure vs. SN explosions, 
iHopkins. Quataert. fc MurravlbOld V Thus, in naive imple- 
mentations of stellar feedback, "star" particles representing 
simple stellar populations (SSPs) distribute "SN energy" 
over neighbouring resolution elements at each time step. 
This procedure is well-known to be inefficient, in that most 
of the thermal energy is radia ted away before it can be con- 
verte d to kinetic energy (e.g. iKatz. Weinberg, fc HernquistI 
Il996l ). 

Three types of sub-grid recipes are commonly 
used to solve the over- cooling problem: injecting the 



energy in kinetic form (e.g. [ Navarro fc Whitd 
Mihos fc HernquistI [Tooi: iKawatal l200ll: iKav et al.1 



Springel fc HernquistI 20031; lOppenheimer fc Dava 



1993 



2002 



2006 



Dalla Vecchia k Schavd l2008l : IPubois fc Tevssierl 12008 
Hopkins. Quataert. Sz Murravl |2012|). s uppressin g radia- 
tive co oling by hand (e.g. iGerritsenI |l997; Mori et al.l 
19971: iThacker fc CouchmanI l2000l: Ikav et all l2002l: 



Sommer-Larsen. Gotz. fc Portinaril l2003l : iBrook et al.l 
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I2OO4I : IStinson et all l2QQ6l : IPiontek fc Steinmetd l201lh . 
and d ecoupling; the di fferen t thermal phases b y hand 



e.g. iMarri fc White! l2003l : IScannapieco et al.1 I2OO6I : 



iMurante et al.l |2010|) . Each solution has its own pros and 
cons and all require the specification of sub-grid parame- 
ters. The different approaches should converge when the 
resolution is increased, although it is not obvious that this 
will in fact happen. 

The inefficiency of thermal feedback is usually at- 
tributed to a lack of resolution: the energy is deposited in 
gas that is too dense, because the hot, low-density, bubbles 
that fill much of the volume of the multiphase ISM are miss- 
i ng;. H owever, as we pointed out in I Dalla Vecchia Sz Schavd 
(|2008l . hereafter DS08), a more fundamental problem is the 
fact that the SN energy is distributed over too much mass, 
which implies that the temperature of the SN heated gas 
is too low and hence the cooling time too short. Indeed, in 
reality one SNII is produced for every ^ 10^ M© of stars, 
and the energy released in the explosion is initially carried 
by ^ 10^ Mq of ejecta. Hence, the ratio of the mass of the 
eject a and the mass of the stellar population that released 
the energy is small <^ 1. In contrast, in simulations this ratio 
is ^ 1, unless the mass of star particles is very large com- 
pared with that of the surrounding gas resolution elements, 
which is typically not the case. Increasing the resolution does 
not change the ratio between the mass of a star particle and 
the mass of the neighbouring resolution elements (although 
it may in the case of grid simulations) and hence is unlikely 
to solve the over-cooling problem by itself. 

The temperature jump of the gas receiving feedback en- 
ergy can be increased by storing the energy until it suffices 
to heat the gas by a desired amount. Indeed, this strategy is 
for example used in some sub-grid recipes for feedback from 
active galactic nuclei (AGN), which store the AGN energy 
in the black hole until the neig;hbouring; g;as can be heated 
to a desired temperature ([Booth Sz Schavd l2009l ). This ap- 
proach is, however, not suitable for SN feedback. Storing the 
energy in a star particle would not help because standard 
implementations of thermal SN feedback are inefficient even 
if all the SN energy of an SSP is released at once. Storing 
energy in gas particles is undesirable because it would make 
the feedback non-local, which would for example mean that 
heavy elements released by the star particle are less likely 
to be carried by outflows. 

An alternative way to increase the temperature jump 
is to reduce the ratio of the heated mass to the mass of a 
star particle. This can be done by reducing the number of 
neighbouring resolution elements that are heated, but that 
does not help if even a single resolution element contains 
too much mass. Moreover, this approach would result in a 
range of temperatures if not all resolution elements have 
the same mass. To guarantee the efficiency of the feedback, 
one would like to be able to specify the temperature jump 
of the gas that receives energy. This can be accomplished 
by making the thermal feedback stochastic: the probability 
that a neighbouring resolution element is heated will then 
depend on the desired temperature jump and on the ratio 
of the mass of the star particle to that of the neighbouring 
gas resolution element. 

A stochastic approach to the rmal feedback was tried, 
and f ound to be effective, by iKav. Thomas, fc Theuna 
(|2003l ) in smoothed particle hydrodynamics (SPH) simula- 



tions of groups of galaxies. For each simulation time- step 
and each star particle, they integrated the energy released 
by SNe and distributed it stochastically to its nearest gas 
neighbour. They only considered the limiting case in which 
the temperature jump was sufficiently high that the number 
of particles that could be heated per time step was less than 
one. 

In this paper w e g;en eralise the method of 
iKav. Thomas, fc Theund (| 20031 ) to work also for tem- 
perature jumps sufficiently low for multiple neighbours 
to be heated. Using SPH simulations of isolated disc 
galaxies embedded in dark haloes with total mass 10^*^ 
and 10^^ h~^ M©, we show that our implementation of 
thermal feedback is able to strongly suppress SF, to alter 
the morphology of the galaxy and to generate galactic 
winds. Reassuringly, for our high-resolution simulations 
we reproduce the results of DS08, where we simulated the 
same disc galaxies with kinetic feedback. 

This paper is organised as follows. We compute the en- 
ergy provided by core collapse SNe (SNII) in Section [21 
where we show that for a standard initial mass function 
(IMF), a single star particle produces enough SN energy to 
heat a gaseous mass equal to the mass of the star particle 
by a few keV, a temperature that is sufficiently high for ra- 
diative cooling to be relatively inefficient. We present our 
numerical implementation for thermal SN feedback in SPH 
in Section O where we also compute some useful quanti- 
ties as a function of the free parameters of the method: the 
amount of SN energy injected per unit stellar mass and the 
desired temperature jump. We also explain how the method 
could be adapted for grid simulations. We dedicate Section |4] 
to the derivation of the resolution constraints of our feed- 
back recipe, which follow from the requirement that the ra- 
diative cooling time exceeds the sound crossing time across 
the heated resolution element so that the thermal energy is 
efficiently converted into kinetic form. After describing our 
simulations of isolated discs galaxies in Section [5l we present 
our results in Section [G] Finally, we summarise and discuss 
our conclusions in Section [71 

Videos and high-resolution images can be found at: 
http : //www . strw . leidenuniv . nl/DS12/ 



2 ENERGY PROVIDED BY SNII 

Each star particle is treated as a simple (or single) stellar 
population. Thus, its stellar content is simply described by 
an IMF, $(M). The number of stars per unit stellar mass 
ending their life as SNII, nsNii, is then the integral of the 
IMF over the mass range [Mo, Mi], 

nsNii= / $(M)dM, (1) 

J Mq 

where Mo and Mi are the minimum and maximum ini- 
tial mass of stars that will explode as core collapse SNe. 
Throughout the paper we will use a Chabrier IMF and 
the mass interval [Mo, Mi] = [6,100] M©, although we 
also report calculations for a Salpeter IMF and for the 
widely used mass range of [8, 100] Mq (6-8 Mq stars ex- 
plode as electr on capture SNe in models with convective 
overshoot; e.g. IChiosi. Bertelli. fc BressanI Il992l ). For the 
Chabrier (Salpeter) IMF we obtain, for the mass range 
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[Mo, Ml] = [6, 100] Mo, nsNii = 1.736 x 10"^ M"^ (1.107 x 
For the mass range [Mo, Mi] = [8, 100] M©, we 



10" 



M-^) 



have nsNii = 1.180 x 10"^ M"^ (0.742 x 10"^ M" 

The total available energy per unit stellar mass provided 
by SNII, esNii = ^snii^snii, is given by 



esNii = 8.73 X 10'^ erg g"' 



^SNII 



1.736 X 10-2 Mq^ 



^51, (2) 



where ^snii = ^51 x 10^"^ erg is the available energy from a 
single SNII event and we will assume E^i = 1. The amount 
of energy from SNII available in a SSP particle is then 
^*esNii, where m* is the initial mass of the star particle. 

If the energy is used to heat a gas mass mg,heat, then 
the corresponding temperature increase is given by 



AT=(7- 



J-j— ; esNii- 



"^gjheat 



.4.23X10^X1 ^^^ r 

\ 1.736 X 10-2 Mq^ 



0.6 



E51 



m* 



"^gjheat 



(3) 



where 7 = 5/3 is the ratio of specific heats for an ideal 
monatomic gas, A;b is the Boltzmann constant, mu is the 
mass of the proton, and /jmu is the mean particle mass. 
We have assumed the gas to be monatomic and neglected 
the energy used to increase the degree of ionisation of the 
plasma. 

The standard SPH approach is to distribute the SN 
energy over all neighbours of a star particle o The heated 
mass is then ?7ig,heat = Nnghrrig, where A^ngb is the number of 
neighbouring particles (typically 32-64, we use 48 in our sim- 
ulations) and ?7ig is the mass of a single gas particle. Assum- 
ing 771* = TTig, we can see from equation (|3]), that the average 
temperature increase for the heated gas particles is '-^ 10^ K, 
which falls in the temperature regjime for which the cool- 
ing tim e is relatively short (e.g. IWiersma. Schave. Sz SmithI 
|2009a|). Note that this procedure leads to even lower tem- 
perature increases if m* < mg, which happens if multiple 
star particles are produced by each star-forming gas parti- 
cle. Heating only a single gas particle would give a temper- 
ature increase of ^ 10^"^ K and, as we will show, a much 
longer cooling time. As we will describe in the next section, 
we will make the temperature increase AT a free parameter 
and stochastically heat neighbouring gas particles. 



3 THERMAL FEEDBACK IMPLEMENTATION 

For simplicity and for consistency with DS08, we assume 
that all the SN energy produced within a star particle be- 
comes available in a single time step. Once a stellar parti- 
cle reaches an age tsN = 3 x 10^ yr, corresponding to the 
maximum lifetime of stars that end their lives as core col- 
lapse SNe, it stochastically injects thermal energy into its 
surroundings. Another reason why we prefer to impose this 
small time delay, is that it prevents the injection of energy 



before the release of heavy elements by (the progenitors of) 
SNe, a process that happens continuously in our stellar evo- 
lution prescription ([Wiersma et al.ll2009bh . 

We note that it is straightforward to modify our imple- 
mentation of thermal feedback to the case where SN energy 
is released stochastically over multiple time steps. All one 
needs to do, is to replace esNii in the equations below by 
the energy per unit mass that is released by the star parti- 
cle over the current time step. We tested this approach and 
did not find any significant difference. 

Before providing a detailed description of our imple- 
mentation of thermal SN feedback, we first describe the two 
free parameters used in our model. 



3.1 Free parameters 

Our recipe uses two free parameters. The first parameter, 
/th, is the fraction of the available SN energy that is actually 
used in performing the feedback [j The value of /th is between 
zero and unity, and can be used to control the efficiency of 
the feedback. Note that this parameter is not particular to 
our method, it is common to all implementations of SNII 
feedback. 

We will present simulations that use either /th = 0.4 or 
/th = 1- The former value is used for comparison with the 
kinetic feedback of DS08, where we u sed the same value. It 
was also used bv ISchave et al.l (|2010|), who already showed 
some results of cosmological simulations that used the pre- 
scription for thermal feedback presented here. For most of 
our runs we will use /th = 1 because it is an interesting 
limiting case and because high values can be justified for 
thermal feedback since we are in fact simulating radiative 
losses (we do not at any point turn off radiative cooling) . 

As discussed in the introduction, we wish to control the 
temperature increase of the heated gas particles in order 
to avoid the regime of short cooling times. We accomplish 
this by making the temperature increase AT, or rather the 
increase of the thermal energy per unit mass Ae, the sec- 
ond free parameter. Although Ae is the parameter we ac- 
tually use, we will often refer to the more intuitive corre- 
sponding value of AT, computed assuming the gas to be 
fully ionised (/i = 0.6). Our fiducial temperature increase is 
AT = 10^"^ K, but we will explore a range of values. 

Simulations employing the recipe described here should 
use values of AT sufficiently high to avoid catastrophic cool- 
ing. As shown analytically in section [4] the required mini- 
mum value will depend on the resolution. The remaining 
freedom, particularly the value of /th, including possible 
dependencies on local physical conditions, can be used to 
calibrate the method by comparing the predictions for the 
quantities of interest to observations. Naturally, the best-fit 
parameter values may depend on other numerical and phys- 
ical parameters. 



^ The standard approach is to weigh the contribution to each re- 
ceiving gas particle by the SPH kernel. We ignore this weighting 
here for simplicity and because differences between SPH neigh- 
bours are by definition poorly resolved. 



^ We do not treat the energy per SN, £^51, as a free parameter 
because it is relatively well constrained by observations and be- 
cause changes in £^51 are in any case degenerate with changes in 
/th- 



© 0000 RAS, MNRAS 000, 000-000 



4 C. Dalla Vecchia & J. Schaye 



3.2 Distributing the energy 

In this section we derive the stochastic formulation for the 
energy injection. The energy released by a single star par- 
ticle is shared among a fraction of the A^ngb neighbouring 
resolution elements. We will hereafter refer to resolution el- 
ements as "particles" but note that the same method will 
also work for grid simulations (replace "particle" by "cell" ) . 
We give each gas particle the same probability p of receiving 
energy from the star, irrespective of its mass (and, for the 
case of SPH, irrespective of its kernel weight). We draw a 
random number < r < 1 for each star-gas particle pair, 
and increase the internal energy of the gas particle by Ae if 
r < p. We will now derive the value of the probability p. 

The expectation value for the total amount of energy 
from SNII injected by a single star particle in the surround- 
ing medium is 



E^- 



^ngb 

: pAe ^ rrii 



(4) 



where Ei and rrii are, respectively, the total energy given to 
and the mass of gas particle i. We require the mean injected 
energy to equal the energy contributed by the star particle, 
fthm^esNU, from which the probability p follows: 



P = /th 



esNii 



m* 



^' E^T"^^ 



(5) 



Thus, the probability that a gas particle is heated is pro- 
portional to /th, the fraction of total available SNII energy 
that the star particle shares with the surrounding gas, and 
inversely proportional to Ae, the amount of thermal energy 
per unit mass that is given to each heated gas particle. 

The stochastic treatment breaks down if the probability 
is larger than one, because in that case the average amount 
of injected energy is lower than the required value. Imposing 
p < 1 puts a constraint on the value of the parameter Ae: 



Ae > /thesNii — Ar~r~ 

2=1 ' 



/th^SNII 



ngb 



(6) 



where the last equality holds exactly if all particles (gas and 
stars) have the same mass, as is usually approximately true 
for SPH simulations. 

For a Chabrier IMF, assuming A^ngb = 48 and /th = 1, 
we obtain Ae > 1.82 x lO^'^ erg g~^, which corresponds to 
a temperature increase AT > 8.8 x 10^ K. For our fiducial 
choice for the temperature increase of 10 "^ K we are well 
above this lower limit. For the case of grid simulations it is 
in principle possible that Ei=f^ '^^^ ^ ^^* ^^^ hence that 
p > 1. In that case it is necessary to increase Ae above the 
chosen value in order to limit the probability to unity. 

The sum of the probability over all neighbours gives the 
expectation value for the number of heated neighbours: 



(A^heat) = fth 



esNii m^Nngh /thCSNii 



Ae 



e:: 



rrii 



Ae 



(7) 



where the last equality holds exactly if all particles have 
the same mass. The average number of heated neighbours 
is inversely proportional to Ae and proportional to /th- Ex- 
pressed in terms of the temperature increase AT, the mean 



number of heated neighbours is 



(iVh, 



1.34 Esi 



nsNU 



1.736 X 10-2 M, 



?(&) 



/ti 



AT 
10''-5 K 



(8) 



Ideally, the energy should on average be shared with at 
least one gas neighbour to make the feedback local to the 
star particle and to ensure that the metals released by mas- 
sive stars can be driven outwards. By injecting all the SNII 
energy from a SSP at once, we can satisfy this constraint 
for temperature increases that are sufficiently large to make 
radiative cooling (initially) inefficient. 



4 ENSURING EFFECTIVE FEEDBACK: 
RESOLUTION REQUIREMENTS 

The thermal feedback can only be effective if the heated gas 
responds hydrodynamically to the temperature increase be- 
fore the thermal energy is radiated away. This implies that 
the sound-crossing time scale across a heated resolution ele- 
ment, ts, must be short compared with the radiative cooling 
time scale in the heated gas, tc. If this condition is satisfied, 
then the gas will start to expand adiabatically, doing work on 
its surroundings and converting thermal energy into kinetic 
energy. The ratio ts/tc can be decreased either by increas- 
ing tc, which can usually be accomplished by increasing the 
temperature, or by decreasing ts, which means increasing 
the temperature and/or the spatial resolution. 

The sound crossing time across a resolution element of 
length h is 



h 



_ h_ _ 1^ IJ^rnu\^ 
"~ cs~ \-fkB J Ti/2 



= 1.15xl0^yr(^) 



1/2 



10^-5 K 



-1/2 



h 



(9) 



100 pc, 

where Cs is the local sound speed. This time scale is related 
to the classical definition of the Courant time step. At = 
Ch/cs = Cts, where C (< 1) is the Courant factor. 
We define the radiative cooling time as 

where u and A are the internal energy and radiative cooling 
rate per unit volume, respectively, and p is the gas density. 
The internal energy per unit volume can be written as 



/cbT 



/cbT 



^ ^ ^-nn, (11) 

7-1 pmu 7-1 M-^H 

where nu and Xu are the hydrogen number density and 
mass fraction, respectively. Before showing results for gen- 
eral cooling functions, we will consider the case of pure 
Brehmsstrahlung, which dominates the cooling rate for 
T > 10^ K and for which the cooling rate is given by 
(IOsterbrocklll989l l: 



A ^ 7.99 X 10 ^^ erg cm ^ s' 



-■(r 



na 



giVeiVmi + ''7HeII + ?7HeIIl) 



: 1.12 X 10"^^ erg cm"^ s"^ (- 



nu 






107.5 K 



107.5 K 



1/2 



1/2 
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Figure 1. Contour plot of the logarithm of the ratio between the 
cooling time and the sound-crossing time over the SPH kernel as 
a function of the density and the temperature of the gas. The 
cooling time scale is computed using tabulated cooling rates for 
primordial (solid contours) and solar (dotted contours) element 
abundances and assuming collisional ionisation equilibrium. We 
show the ratio for the particle mass used in our simulations of the 
10-*^^ h~^ Mq halo. The vertical dashed line marks the fiducial 
heating temperature of 10^'^ K and the horizontal dashed line 
the star formation density threshold in our simulations. 



(1 + Xh)(1+3Xh) 
8X2 



(12) 



where 771 = m/nu is the number density of species i relative 
to hydrogen and we assumed the plasma to be fully ionised 
in the last step. We assumed for simplicity that the Gaunt 
factor Qi pi 1.4 (consist e nt wi th the approximated equation 
given by iTheuns et al.l (| 19981 ) for our fiducial temperature 
increase) and a primordial composition. Hence, we obtain a 
radiative cooling time of 



'. 3.26 X 10'^ yr 



( ^^ rV 

vicm-37 V 



107.5 K 



1/2 



VO.6/ V 0-13 J ' 



where 



Xh(1 + Xh)"'(1 + 3Xh)"' 



(13) 



(14) 



and /(Xh = 0.752) ^ 0.13. 

Thus, we can write for the ratio between the two time 
scales 



- = 2.; 

ts 



X 10^ 



VI cm-3; 



lO^-s Kj VlOO pc 



h 



/ /i X-3/2 //(X. 

V0.6/ V 0.1; 



V 0.13 



(15) 



For the case of SPH simulations, the spatial resolution 



is given by the gas particle's smoothing kernel which can be 
approximated as 



3 e:=t 



ngb ^ 



47r 



1/3 



3 A^ngb (m) 
4:7T munu 



Xu 



1/3 



(16) 



where (m) is the average mass of the gas particles. Substi- 
tuting the above approximation into equation ([15)) we obtain 
the ratio of time scales 



^-98(-^) 
is VI cm -^ / 



N, 



ngb 



48 



-1/3 



-2/3 



vo.6y 



(m) 



10^-5 Kj V7x 10^ Mo 

3/2 r gjXu] 
0.14 



-1/3 



(17) 



X^^^'^ /{Xh) and we used the particle mass 



where g{Xu) — X^ 

appropriate for our simulations of the 10^^ h~^ Mq halo. 

We expect cooling losses to be important for tc ^ ts. 
The exact value of the ratio ft = tc/ts required to ensure 
efficient feedback can only be determined using simulations, 
but we expect it to be similar to 10 and will therefore use 
this as our fiducial v alue. This agre e s well with the recent, 
independent work of ICreasev et al.l (|201l|), who derived a 
resolution criterion for shock capturing in SPH and adap- 
tive mesh refinement (AMR) simulations. Their criterion is 
based on comparing the rates of shock heating and radia- 
tive cooling in a shock front, and ensures that shock heating 
overwhelms cooling in order to avoid numerical over-cooling. 
Following the suggestion of DS08, they applied their results 
to the injection of thermal energy, and obtained a criterion 
that basically translates into tc/ts > 8lj 

From the relation between the sound crossing time and 
the Courant time step, one can estimate the number of simu- 
lation time steps it would take the gas to radiate its thermal 
energy if it cooled isochorically: ristep ^ tc/At = ftts/At — 
ft/C (^ 30 for the fiducial value of ft and C = 0.3). In most 
implementations of SPH the sound speed in t he Courant 
criterion is replaced by a "signal velocity" (e.g. iMonaghanI 
1997), t'sig ^ 2cs, which would more than double the number 
of time steps. 

Inverting equation (pT|) , we find the following maximum 
density for which the feedback is expected to be effective. 



nu,tc^ftts = 31 cm 



3/2 / n \ -3/2 
Jt 



107.5 K 

/ \ \ -1/2 



7 X 10^ M 



0.67 I 0.14 



iVngb 

48 

3/2 



-1/2 



(18) 



The critical density is thus proportional to T^/'^ {m)~ ' . 

For the case of AMR simulations, the spatial resolution 
is given by the linear size of the grid cell, which is commonly 
decreased until it is some factor, /j, smaller than the local 
Jeans length: h < Lj/fj. Expressing the Jeans length as a 
function of temperature and density, L j = Cs ^7v/{Gp), and 



^ Note that ICreasev et al.l (|201lh derived their criterion by esti- 
mating the velocity of a blast wave at the blast radius of half the 
mean inter-particle distance. 
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substituting this into equation ([15} , we obtain 



>50 



Vlcm-3y VlO^-^K 

V0.6/ V 0.15 



To 



104 K 



-1/2 



/j 



(19) 



where To is the initial gas temperature (from which the Jeans 
mass, hence the spatial resolution, is derived), ^/(Xh) = 
X^ ' /(Xh), and we assumed that the Jeans length is re- 
solved with four resolution elements. The equality holds for 
h = Lj/fj. Note that the same relation applies to SPH sim- 
ulations with particle mass at least /| times smaller than 
the Jeans mass in gas with density nu and temperature To. 
Inverting equation ([T9|), we obtain 



nH,tc=/tts > 25 cm" 



(JLY 
vo.6y 



10^-5 K 
0.15 



104 K 



-1 / n \ 2 
/j 



(20) 



which is the analogue of equation ([18} for AMR (and for 
SPH simulations with particle mass at least /| times smaller 
than the Jeans mass in gas with density nu and temperature 
To). 

We show in Fig.[T]a contour plot of the time-scale ratio 
tc/ts in the temperature-density plane for the fiducial gas 
particle mass rrig = 5.1 x 10^ h~^ M©. We computed tc us- 
ing!: the tabulated values of the collisional cooling rates of 
(jWiersma. Schave. fc Smith 2009a '). which are also used in 
the SPH simulations described below. We combine in the 
same plot the two cases of primordial (solid contours) and 
solar (dotted contours) chemical compositions. The verti- 
cal dashed line marks the fiducial heating temperature of 
10^"^ K, while the horizontal dashed line indicates the den- 
sity threshold for SF in our simulations. For T > 10^"^ K the 
difference between primordial and solar metallicity is very 
small, but at lower temperatures metals reduce the ratio 
tc/ts. Once the temperature has dropped to 10^ K for the 
case of solar abundances (or to values smaller than 10^ K 
for primordial abundances), collisional excitation processes 
cause the cooling rate to increase as the temperature drops 
(e.g. [Wiersma, Schaye, & Smith 2009a). For such tempera- 
tures the equations above, which assumed the cooling rate to 
be dominated by Brehmsstrahlung, will underestimate the 
radiative losses and will therefore overestimate the minimum 
density for which the feedback is expected to be efficient. 

Interestingly, for purely adiabatic expansion, i.e. p oc 
T"^/^, the ratio tc/ts remains constant for the case of SPH 
(eq. [13 )• In that case the gas will follow a track paral- 
lel to the high temperature part of the contours in Fig. [T] 
and tc ex: p~'^'^ . Hence, if radiative losses were unimportant 
initially, so that the hot bubble will start to expand adiabat- 
ically, then radiative losses will remain unimportant as long 
as the cooling rate is dominated by Brehmsstrahlung. For 
AMR, on the other hand, the ratio tc/ts oc T^l^ during the 
adiabatic phase, which implies that radiative losses may al- 
ready become important while Brehmsstrahlung dominates. 



5 SIMULATIONS 

We ran simulations of isolated disc galaxies embedded 
in dark matter haloes with total masses of 10^° and 
VS^^hr^ Mq, where h — 0.73. T he initial conditions are as in 
ISchave h Dalla Vecchial (|2008l ) and DS08, thus the models 
do not include gaseous haloes and all the gas is initially in 
the discs. We adopted a larger gravitational softening length 
than in the previous works for the massive galaxy. We also 
modified the original code as described below. 



5.1 Code and initial conditions 

We use a modified version of the TreePM/SPH code gad- 
get (|SDringelll2005l ) for all the simulations presented in this 
paper. 



We employ the SF recipe of I Schave h Dalla Vecchial 
(|2008l ). which we briefiy describe here. Gas denser than the 
critical density for the onset of the thermo-gravitational in- 
stability (riH ^ 10~^ — 10~^ cm~^) i s expected to be multi- 
phase and star- forming (Schaye 2004). We model such gas by 
imposing a minimum temperature fioor given by an effective 
equation of state with pressure P oc p^g""^^ for densities ex- 
ceeding riH = 0.1 cm~^, normalised to P/k = 10^ cm~^ K at 
the threshold. We use 7eff = 4/3 for which both the Jeans 
mass and the ratio of the Jeans length and the SPH ker- 
nel are independent of the density, thus preventing spurious 
fragmentation due to a lack of numerical resolution. 

We introduce a different definition of star-forming par- 
ticle. In the previous works a gas particle was fiagged as star- 
forming if it crossed the density threshold n^ = 0.1 cm~^ 
while its temperature was below T = 10^ K. The particle 
then remained star-forming until its density fell below the 
threshold density or the particle was promoted to a wind 
particle. 

In the present work we proceed as follows. Each parti- 
cle is free to cool to lower temperatures, but not below the 
temperature TeoS imposed by the effective equation of state. 
A gas particle is star-forming if 



^ logio T < logio ^Eos + A logio ^EoS 



(21) 



where A log^^Q TeoS is a free parameter which we set to 
AlogioTEoS = 0.5 dex. 

The Kennicutt-Schmidt SF law is analytically converted 
an d implemented as a pre s sure l aw. As we demonstrated 
in I Schave fc Dalla Vecch ia (2008), our method allows us 
to reproduce arbitrary input SF laws for any equation of 
state with o ut tu ning any parameters. We use the observed 
iKennicutti (Il998l ) law 



E* - 1.5 X 10"^ M0yr"^kpc"^ 



1.4 



. 1 M0 pc 



-2 



(22) 



where the different normalisation accounts for the fact that 
we are using a Chabrier IMF. 

Radiative cooling and heating were included 
using; tables for hyd r ogen and helium from 
I Wiersma. Schave, fc SmithI (|2009al ) . The cooling tables 
were generated using the pub licly available package CLOUDY 
(version 07.02; iFerlandl l200p|) . assuming ionisation equilib- 
rium in the presence of the iHaardt fc Madaul (|200lh model 
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Table 1. Simulation parameters: total mass, M^aio; fraction of SN energy injected, /th; temperature jump, log^o AT; total number 
of particles, A^tot; total number of gas particles in the disc, A^discJ mass of baryonic particles, mb; mass of dark matter particles, 
^^dm; gravitational softening of baryonic particles, Cb; gravitational softening of dark matter particles, cdm; thermal feedback included, 
(Feedback). Values different from the fiducial ones are shown in bold. 



Simulation 


Mhalo 
(/i-l M0) 


/th 


logio AT 
(K) 


A^tot 


iVdisc 


(^-1 M0) 




(h-^ pc) 


(h-^ pc) 


Feedback 


GIO-NOFB 


IQlO 


— 


— 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


10 


17 


N 


G 10-040-70 


IQlO 


0.4 


7.0 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


10 


17 


Y 



GlO-100-65 
GlO-100-70 
G 10- 100-75 
GlO-100-80 
GlO-100-85 



1.0 


6.5 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10"^ 


10 


17 


Y 


1.0 


7.0 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


10 


17 


Y 


1.0 


7.5 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


10 


17 


Y 


1.0 


8.0 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


10 


17 


Y 


1.0 


8.5 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


10 


17 


Y 



G10-100-75-LR08 
G10-100-75-LR64 


10^0 
10^0 


1.0 
1.0 


7.5 
7.5 


625 061 
78 132 


29411 
3 676 


4.1 X 10^ 
3.3 X 10^ 


1.9 X 10^ 
1.5 X 10^ 


20 
40 


34 
68 


Y 
Y 


G12-NOFB 


10^2 


— 


— 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


N 


G 12-040-70 


10^2 


0.4 


7.0 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


Y 



G12-100-65 


10^2 


1.0 


6.5 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


Y 


G 12- 100-70 


10^2 


1.0 


7.0 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


Y 


G 12- 100-7 5 


IOI2 


1.0 


7.5 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


Y 


G12-100-80 


10^2 


1.0 


8.0 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


Y 


G12-100-85 


IOI2 


1.0 


8.5 


5 000 494 


235 294 


5.1 X 10^ 


2.4 X 10^ 


46 


79 


Y 



G12-100-75-LR08 
G12-100-75-LR64 



10^ 
10^ 



1.0 
1.0 



7.5 

7.5 



625 061 
78 132 



29411 
3 676 



4.1 X lO"* 
3.3 X 10® 



1.9 X 10® 
1.5 X 10''' 



92 

184 



158 
316 



Y 

Y 



for the 2; = UV background radiation from quasars and 
galaxies, and the cosmic microwave background. 

We implement ed a version of the tira e -stepp ing algo- 
rithm described bv lPurier k, Dalla Vecchial ((2013). Inactive 
particles that receive feedback energy are immediately acti- 
vated so that they can respond promptly to their new en- 
ergetic state. T heir and their ac tive neighbours' signal ve- 
locities (see e.g. lMonaghanlll997l ) are also updated in order 
to calculate the size of the next time-step consistently. The 
new time-step is then pr opagated to the inactive neighbours 
following the scheme of Saitoh fc M akino (2009). We refer 
the reader to iDurier fc Dalla Vecchial (|2012l ) for a detailed 
discussion of the benefits introduced by the scheme. 

However, we do not expect the integration accuracy of 
our simulations to i mprove significantly compared with e.g. 
DS08. As argued bv lPurier fc Dalla Vecchial (|2012l ). impos- 
ing a limit to the maximum allowed time-step (e.g., by signif- 
icantly reducing the gravitational softening, as we have done 
here and in DS08), may also maintain good energy conser- 
vation in the case of strong energy perturbations. Indeed, we 
found only small deviations in the global properties of the 
outflows when running the same fiducial models without the 
time-step limiter. 

The initial conditions are based on the model of 
ISpringeh Di Matteo, fc HernquistI (|2005l ) and are described 
in DS08. The model consists of a dark matter halo, a stellar 
bulge, and an exponential disc of stars and gas. The circular 
velocities at the virial radii are 35.1 and 163 kms~^ for the 
10^° and 10^^ /i~^ M© haloes, respectively. The virial radii 
are 35.1 and 163 h~^ kpc. The halo is rotating and has a 
dimensionless spin parameter A = 0.33. The disc contains 4 
percent of both the total mass and the total angular momen- 
tum. The bulge contains 1.4 percent of the total mass and 
has a scale length one tenth of that of the disc. The bulge 



has no net rotation. The initial gas fraction of the disc is 30 
percent, the remaining 70 percent of the disc mass is made 
up of stars. The vertical distribution of the stellar disc has 
a constant scale height of 10 percent of the radial disc scale 
length. 

Except for our low-resolution runs, the total number of 
particles in each simulation is 5,000,494, of which 235,294 are 
gas particles in the disc. The baryonic particle mass for the 
10^°/i"^ Mo (lO^^/i"^ Mo) halo is mb = 5.1 x lO^/i"^ Mo 
(mb = 5.1 X 10"^ /i~"^ Mo)- The gravitational softening 
length was set to eh — 10/i~^ pc for the baryons and to 
(mDM/mb)"^^^eb ~ I7h~^ pc for the dark matter in the 
10^^ h~^ Mo halo. The softening for the massive galaxy is 
scaled up in proportion to rrij^^. 



5.2 Simulation parameters 

Table [1] lists the simulations we have performed. Each sim- 
ulation was evolved for 500 Myr. The simulations are la- 
belled with the prefix GIO and G12 for the 10^° and the 
10^^ h~^ Mo haloes, respectively, followed by the percent- 
age of the SN energy that is injected and by the logarithm 
of the temperature increase. For example, GlO-040-75 refers 
to the 10^° /i~^ Mo halo, with the SN feedback injecting 40% 
of the total available energy and increasing the gas particle 
temperature by AT = 10^"^ K. We have run several vari- 
ations of the fiducial models, GlO-100-75 and G12-100-75. 
The list of all the models follows: 

• One run without SN feedback {G[10,12]-NOFB). 

• One run injecting 40 percent of the available SN en- 
ergy, G[10,12]-040-75. This is used for comparison with the 
kinetic feedback simulations of DS08, which also used 40 
percent of the energy. 
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AT = lo^-^ K 



AT = 10^-5 K 



AT = lO^-^ K 




Figure 2. Projections of the gas density (top row: edge-on; middle row: face-on) and temperature (bottom-row: edge-on) for models 
GlO-100-65 (left column), the fiducial model GlO-100-75 (middle column) and GlO-100-85 (right column) at time t = 250 Myr. The 
white arrows in the top row show the velocity field. Images are 17.5h~^ kpc on a side. The colour coding is logarithmic in density 
(—4.3 < log^o ^H/cm""^ < —0.3) and temperature (3.7 < log^o ^/K < 5.8). The colour scale is indicated by the colour bars in each 
column. 



• One set of runs varying the temperature increase 
logioAT = [6.5,7.0,7.5,8.0,8.5]. These runs are labelled 
G[l 0,1 2] -100- [6 5, 70, 75, 80, 85], respectively. 

• Two runs in which the number of particles was de- 
creased by factors of 8 and 64, respectively {G[10,12]-100- 
75-LR08, G[10,12]-100-75-LR64)- 



6 SIMULATION RESULTS 

In this section we describe tests of our implementation of 
thermal feedback. We first show the effect of varying the 
temperature increase. We proceed with a comparison to the 
kinetic feedback model of DS08. We conclude the section by 
reporting the results of resolutions tests. 



To eliminate differences other than the feedback recipe 
from the comparison with DS08, we repeated the simulations 
of models ml2 and ml2nowind of DS08 with the new version 
of the code. We changed the softening to the one used in this 
work, and used the time-step limiter for the feedback run. 



6.1 Dependence on the temperature increase 

We ran a set of simulations with /th = 1 and varying 
logio ^T = [6.5, 7.0, 7.5, 8.0, 8.5] to study the dependence 
on the temperature increase. We first describe the morphol- 



© 0000 RAS, MNRAS 000, 000-000 



Galactic outflows with thermal SN feedback 9 



AT = lO^-^ K 



AT = 10^-^ K 



AT = lO^-^ K 




Figure 3. Projections of the gas density (top row: edge-on; middle row: face-on) and temperature (bottom-row: edge-on) for models 
G12-100-65 (left column), the fiducial model G12-100-75 (middle column) and G12-100-85 (right column) at time t = 250 Myr. The 
white arrows in the top row show the projected velocity field. Images are 4:5h~^ kpc on a side. The colour coding is logarithmic in 
density (—4.3 < log^o^H/cm"^ < 1.3) and temperature (4.3 < log^o^/K < 6.3). The colour scale is indicated by the colour bars in 
each column. 



ogy of the galaxies and their outflows, and proceed with 
discussing the SF histories and the outflow properties quan- 
titatively. 



6. 1.1 Morphology 

For each halo, we compare here the fiducial model, which 
has log^o ^^ — '^•^5 ^o our most extreme models, which use 
logj^Q AT = [6.5, 8.5]. Fig.[2]shows projections of the density 
(top and middle rows) and temperature (bottom row) for the 
dwarf galaxy at time t = 250 Myr for three different values 
of the temperature increase (from left to right, log^g ^^ — 



[6.5, 7.5, 8.5]). Recall that all models inject the same amount 
of energy per unit stellar mass. 

The morphology of the galaxy is irregular in all cases. 
As the heating temperature is increased, the low-density 
bubbles in the ISM and in the circumgalactic medium in- 
crease in size and open up vertical channels through which 
the outflows can escape. Consequently, the outflow becomes 
more collimated along the vertical axis as AT increases, en- 
hancing its bipolarity. The velocity field overplotted in the 
top row shows that the flow is faster for larger AT, while 
the bottom row shows that the outflowing gas is also hotter. 

Fig. [3] shows edge-on projections of the density (top 
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AT = 10^-^ K 



AT = lO"^-^ K 



AT = 10^-^ K 
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-4-2 2 

logio iIh [cm"^] 



-4-2 2 

logio iIh [cm~^] 



Figure 4. Two-dimensional mass-weighted probability density distribution of the gas within 0.2rvir of the dwarf galaxy in temperature- 
density space at time t = 250 Myr. The colour coding indicates /g = (dM/M)/d log^^Q nn/dlog^o T. The three panels correspond to the 
same three models as were shown in Fig. [2] The twelve contours of /g (black, dotted) are equally spaced in the range showed by the colour 
bar. The red, dashed curves indicate radiative cooling time contours and their labels indicate log2o(^c/yr)- Below the solid, green line 
photo-heating by the UV background dominates over radiative cooling. The horizontal, dashed line indicates the heating temperature 
AT. The vertical, dotted line marks the threshold for star formation. The imposed, effective equation of state is visible for densities 
larger than the threshold density. 
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Figure 5. As Fig. ID but for the models of the massive galaxy shown in Fig.O 



row) and temperature (bottom row) for the massive galaxy 
for the same three different values of the temperature in- 
crease. The dependence on AT is more evident for this more 
massive galaxy. For AT — 10^"^ K (left column), most of the 
outflowing gas is confined to a region around the disc. The 
disc looks puffed up as the expelled gas is deposited just 
outside it, and the gas breaks up in little blobs. The veloc- 
ities are small and the velocity field does not show a clear 
preferential orientation. There is a galactic fountain, but no 
large-scale galactic wind. 



In contrast, the fiducial model (AT = 10^"^ K; middle 
column) shows a clear bipolar outflow, which is sustained 
until the end of the run. The outflow is mostly driven from 
the inner part of the disk where a large fraction of the SF is 
taking place, and it is collimated by the disc which impedes 
motion within its plane. The temperature map shows several 
cold blobs above and below the galactic disc. The blobs come 
from the disc and are moving outward, thus the outflow is 
ejecting parcels of cold gas. Cold blobs are seen falling back 
onto the disc at large radii. 
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Figure 6. The Kennicutt-Schmidt SF relation for a selection of 
GIO (left panel) and G12 (right panel) models at t = 250 Myr. 
Surface densities are computed in cylindrical annuli containing 
a constant number of particles and including all particles with 
vertical coordinate |2;| < 2 kpc. The tilted line shows the observed 
Kennicutt SF law, equation (|22|) , whereas the vertical line shows 
the gas surface density below which SF is observed to bec ome 
inefficient (for more details see lSchave &: Dalla Vecchiall2008l V All 
models are in excellent agreement with the observations. 



The highest AT run (AT = 10^"^ K; right column) 

looks similar to the fiducial model, but the wind is hotter 
and moving faster and the circumgalactic medium contains 
fewer cold blobs. 

Videos illustrating the time evolution of the 

model galaxies are available at this web address: 
http : //www . strw . leidenuniv . nl/DS12/ 



6.1.2 Gas phase distribution 

The mass-weighted probability distribution function (PDF) 
of the gas at t = 250 Myr in the (nn , T) plane is shown in 
Figs. [Hand [5] for the models shown in Figs. [2] and [S] respec- 
tively. To limit the effects of the vacuum boundary condi- 
tions (mainly adiabatic cooling to extremely low tempera- 
tures and densities), we include only the gas inside a spher- 
ical volume of radius 0.2rvir, but we normalise the PDFs 
to the total gas mass. The vertical, dotted line marks our 
threshold density for SF, and at higher densities the im- 
posed, effective equation of state is clearly visible. Contours 
of constant radiative cooling time (dashed, red lines) are 
over-plotted 1j The contour labels indicate integer values of 
log^Q(tc/yr). The equilibrium between radiative cooling and 
photo- heating by the UV background is shown by the solid, 
green curve. Below this curve radiative heating dominates 



^ Note that the tc contours have been calculated assuming a con- 
stant mean molecular weight, ji = 0.6. 



over radiative cooling. Note that gas can reach tempera- 
tures lower than the equilibrium value if it cools adiabati- 
cally. In the absence of feedback, none of the gas below the 
SF threshold would have temperatures significantly above 
the equilibrium value. Finally, the horizontal, dashed (green) 
lines indicate the feedback heating temperatures AT. Gas 
in which feedback energy has been injected initially resides 
near these lines and, as long as radiative cooling is unim- 
portant, will expands adiabatically, exiting the star- forming 
region at different temperatures. Indeed, the phase diagrams 
show different distributions of shock-heated gas for different 
AT, with the low-density, high-temperature regions being 
more populated for increasing AT. 

For the dwarf galaxy (Fig. [J) the temperature-density 
distributions are similar for all values of AT, although there 
are some differences in the high temperature regime. We 
will show later that the SF histories and the outflows are 
also very similar for all values of AT and that, given the 
resolution, this is in accord with the results of section [H 
There is very little hot gas. The surface density of the disc, 
and hence the pressure, is too low to confine the heated gas, 
which therefore immediately blows out of the disc and cools 
adiabatically. 

For the massive galaxy changing AT has a much more 
dramatic impact on the distribution of shock-heated gas 
(T > 10^ K). For the model with the lowest value of AT 
(left panel of Fig. [5|), the peak of the PDF lies within the 
density range 10""^ to 10~^ cm""^. This confirms the qual- 
itative result shown in Fig. O most of the outflowing gas 
accumulates in a region around the disc. If we increase AT 
(middle and right panels), the peak in the PDF moves to 
lower densities (^ 10~^ cm""^) thanks to the development 
of a large-scale outflow that moves gas away from the disc. 
The total fraction of shock- heated gas decreases with the AT 
because less gas resides within 0.2rvir if the wind velocity is 
higher. 



6.1.3 Star formation history 

Before discussing the SF histories, we show the predicted 
Kennicutt-Schmidt SF relations in Fig. \6\ The left (right) 
panel shows the same three different models as were shown 
in Fig. [2] (Fig. [3|) for the dwarf (massive) galaxy. Gas mass 
and SF surface densities were computed in annuli contain- 
ing a constant number of gas particles and including all par- 
ticles with vertical coordinate \z\ < 2 kpc. The observed 
Kennicutt-Schmidt law (eq. (|22p : tilted line) and the steep- 
ening at the SF threshold density (vertical line) are well 
matched. This success is not un e xpect ed, as we already 
showed in ISchave fc Dalla Vecchial (|2008l l that the observed 
SF law can be implemented directly in the form of a pres- 
sure law and that this enables the simulations to reproduce 
the observations without the need to tune any parameters 
and irrespective of whether strong feedback is present. While 
feedback determines the surface density of the gas, it does 
not affect the efficiency of star formation at a fixed surface 
density in our models, since pressure and surface density are 
closely related in self- gravitating systems. 

The dependence of the SF histories on AT is shown in 
Fig. [71 For the dwarf galaxy (left panel), the star formation 
rate (SFR) drops sharply within the first 100 Myr due to 
the strong feedback produced by the initial burst of SF, and 
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Figure 7. SFR as a function of time for the galaxies in the 10"'^*-' and the 10^'^ h~^ Mq haloes (left- and right panels, respectively). We 
vary the temperature increase due to feedback events in steps of 0.5 dex over the range log^Q AT/K = [6.5, 8.5]. The dwarf galaxy's SF 
history is insensitive to the value of AT, while for the massive galaxy the feedback is less efficient for the lowest heating temperature. 



remains nearly constant thereafter. The fact that the sharp 
drop is due to feedback can be seen by comparing to the 
nearly horizontal black, dotted curve, which shows the SF 
history for a run without feedback. The factor by which the 
feedback reduces the SFR is insensitive to AT. This is in ac- 
cord with the calculations presented in Section [H The dwarf 
galaxy has a low surface density and forms most of its stars 
at densities close to the SF threshold of nn = 10~^ cm~^ 
(see Fig. \^ . For such densities and primordial abundances 
we expect cooling losses to be small even for heating tem- 
peratures as low at 10^"^ K. Indeed, for our particle mass of 
7 X 10^ M0 and a heating temperature of T = 10^"^ K, equa- 
tion (|18p tells us that cooling losses should be unimportant 
for densities nn < 10 cm~^. 

However, as Fig. [T] shows, the situation could be differ- 
ent for solar abundances. While the cooling time is insen- 
sitive to the metallicity for T>10^-^ K, for T = 10^-^ K 
it is about an order of magnitude smaller for solar metal- 
licity than it is for primordial abundances. A factor ten in- 
crease in the cooling rate would reduce the maximum den- 
sity for which the feedback is expected to be effective to 
riH ~ 0.3 cm~^ (eqs. [17] and [H]). We would therefore ex- 
pect a significant fraction of the feedback energy to be radi- 
ated away before it can be converted into kinetic form, if we 
were to run the dwarf galaxy simulation with AT = 10^'^ K 
and solar metallicity. Indeed, we have performed such a run 
(not shown) and find the SFR to be much higher. After 
400 Myr it is about 0.04 M© yr~^ which is closer to the run 
without feedback than to the run with the same AT but 
primordial abundances. 

For the massive galaxy the decline in the SFR is more 
gradual (right panel). After a few hundred Myr, all runs 
predict roughly the same SFR except for the one adopting 
AT = 10^-^ K, which has a substantially higher SFR. The 
morphological comparison (Fig. [3|) shows that in this run 
the gas is unable to escape to large radii. Instead, it ac- 
cumulates around the disc, and eventually falls back onto 
it. This is expected, because equation ([T8|) shows that ra- 
diative losses will become important at densities that are 
10 times lower than for the dwarf galaxy, because the par- 
ticle mass is 100 times higher. Moreover, the densities in 
the ISM are higher in the massive galaxy, with many star 
particles forming in gas with densities nn '^ 1 — 10 cm""^ 



(Fig. O. Comparing these numbers to the maximum den- 
sity for which cooling losses are small at this resolution, 
nu < 31 cm-3 (T/IO'^-^)^/^ (eq. [H]), we expect small 
changes for AT = 10^ K and a substantial reduction of 
the feedback efficiency for AT = 10^"^ K. 

Models G12-100-70 and G12-100-75 have similar SF 
histories, suggesting convergence of the results. However, 
models G 12- 100-80 and G 12- 100-85 have SFRs that are 
larger than that of the fiducial model by factors of :^ 25 and 
40 percent, respectively. The trend for the largest AT's likely 
arises from poor sampling of the distribution of SN energy 
in the disc. Indeed, the expectation value for the number of 
heated neighbours decreases with AT, and is (from eq. [7]) 
0.42 and 0.13 for G 12- 100-80 and G 12- 100-85, respectively. 
This shows the importance of locally linking the feedback 
events with the formation of star particles by depositing the 
star particle SN energy into at least one of its neighbours. 
The effect is less severe for the dwarf galaxy because star 
formation is restricted to smaller scales, these scales are re- 
solved with many more particles, and it is easier to eject gas 
in this case due to the lower ISM pressure and the shallower 
gravitational potential well. 

Finally, we note that the SFRs are lower than for the 
kinetic feedback runs of DS08. This difference is, however, 
due to the different amount of energy injected rather than to 
the manner in which this is done (i.e. thermal vs. kinetic). 
We inject more energy here (/th = 1 whereas DS08 used 
/th = 0.4), thus a stronger quenching is expected. In Sec- 
tion 16.21 we will show that the two methods are in fact in 
good agreement with each other. 



6.1.4 Mass outflow rate and wind velocity 

In this section we measure the wind velocity and mass load- 
ing as a function of time and radius. We first briefiy describe 
the method, which is identical to the one employed in DS08. 
We discretise the following integral equation of the net 
mass outfiow rate through a surface S: 



M 



I pv • dS, 

J s 



(23) 



where p is the gas density and v is the gas velocity at any 
position on the surface. Given a spherical shell of radius r 
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Figure 8. Mass outflow rate (left column) and average outflow velocity (right column) measured through a spherical shell at radius 
r = 0.2rvir as a function of time (top row) and at t = 500 Myr as a function of radius (bottom row) for models G10-100-[65,70,75,80,85]. 
The dotted line in the top-left panel indicates the SFR of model GlO-100-75. All other curves are labelled in the legends. 



and thickness Ar centred on the origin, the above equation 
becomes 



M{r, Ar) 



-, -'^'shell 

— y 

Ar ^^ 



. / rriiVi • — , 
Ar ^ n 



(24) 



where A^sheii is the total number of particles within the shell, 
and rrii and r^ are their mass and position, respectively. We 
consider all particles and use Ar = rvir/150. 

The average outflow velocity is also given in discrete 
form, and is the mass- weighted, average radial velocity: 



(v) (r, Ar) 



E 



N.hell , 



E-'^'shell 
2=1 



(25) 



where v^ is the particle velocity. We consider only particles 
moving away from the origin, as indicated by the subscript 
"+" in the above equation. 

To investigate the dependence of the outflow properties 
on the temperature increase AT, we show in Figs. [8] (dwarf 
galaxy) and [9] (massive galaxy) the mass outflow rate (left 
column) and the average outflow velocity (right column). 
The top row shows the evolution of the wind at radius r = 
0.2rvir and the bottom row shows the dependence on radius 
at time t = 500 Myr. 

We first consider the dwarf galaxy (Fig. [8]). Apart from 
the peak velocity, which is reached after only a few tens of 
Myr, all models look very similar. Apparently, the proper- 
ties of the wind are nearly independent of the temperature 
increase AT. As discussed in the previous section, this is 
expected for such a high resolution and for a primordial 



composition. However, for solar abundances the equations 
derived in Section |4] predict that radiative losses do become 
important for AT = 10^"^ K. Indeed, we find that for this 
heating temperature the mass outflow rate is reduced by 
more than an order of magnitude if we assume the metallic- 
ity to be solar (not shown). 

The wind is highly mass- loaded, with mass outflow rates 
at 0.2rvir that are a factor of 10 — 100 higher than the SFR 
(the SFR of the fiducial model is shown as the black, dotted 
curve in the top- left panel). The mass fiux builds up quickly 
in the first 100 Myr and declines gradually thereafter. 

The top-right panel shows that the wind velocity peaks 
at different values at the beginning of the simulations, but 
converges to similar values after ^ 50 Myr. The peak outflow 
velocity depends strongly on the temperature increase, with 
the smallest (largest) value of '^ 400 kms""*^ (^ 600 kms""*^) 
corresponding to the smallest (largest) AT. The peak ve- 
locity is a measure of the average velocity of gas particles 
that are able to move freely to large radii. Because our sim- 
ulations start without a gaseous halo, it is not clear how 
meaningful the early evolution is, given that the high wind 
velocities would have resulted in strong shocks if a gaseous 
halo had been present. However, because the winds flll the 
haloes with gas, the results quickly become insensitive to the 
artiflcial initial conditions. After about 100 Myr the wind 
velocities of all models converge at about 100 kms~^ after 
which they decline gently to several tens of kilometres per 
seconds. 

The convergence in the SFR, mass outflow rate and 
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Figure 9. Mass outflow rate (left column) and average outflow velocity (right column) measured through a spherical shell at radius 
r = 0.2rvir as a function of time (top row) and at t = 500 Myr as a function of radius (bottom row) for models G12-100-[65,70,75,80,85]. 
The dotted line in the top-left panel indicates the SFR of model G12-100-75. All other curves are labelled in the legends. 



outflow velocity for different values of AT indicates that 
similar amounts of star- forming gas has been extracted from 
the disc over time. The plots in the bottom row of Fig. [8] 
confirm that. At each radius the outflow rates and velocities 
are similar for all models, and the gas can efficiently escape 
into the halo, and eventually beyond the virial radius. 

As was already noted by DS08, the fact that the wind 
velocity becomes proportional to the radius as we move away 
from the galaxy (bottom-right panel) is due to travel time 
effects. Because the wind has only been blowing for a finite 
amount of time t, only gas with a velocity greater than 



V — 9% kms 



10 kpcy Vio^ yr 



t 



(26) 



would be present at radius r if the wind velocity is constant 
with radius. Although the wind may in reality accelerate or 
decelerate, our results suggest that travel time effects may 
well be important too. Hence, wind velocities that are ob- 
served to increase with distance do not necessarily imply 
that the wind is accelerating. 

Fig.[9]shows that the picture is different for the massive 
galaxy. The initial wind velocity is sensitive to the heating 
temperature and determines the time at which the outflow 
first passes through the shell at r = 0.2rvir (corresponding to 
the sharp rise in the top panels). For AT > 10^"^ K the mass 
outflow rates converge at about five times the SFR (the SFR 
of model G12-100-75 is indicated by the dotted curve in the 
top-left panel), so the winds are much less mass- loaded than 
for the dwarf galaxy. For lower heating temperatures the 



mass flux is substantially lower. The outflow rate of model 
G12-100-65 even becomes negative after about 320 Myr, in- 
dicating net infall (dotted curve in the top- left panel). The 
wind cannot escape the inner region of the halo, and is con- 
fined within a fraction of the virial radius. This model there- 
fore also predicts a much higher SFR (Fig.fT]). The fact that 
the mass flux drops strongly for AT < 10^"^ K can be under- 
stood by noting that the gas in the central regions and spiral 
arms has densities nn '^ 1 — 10 cm""^ (Fig-E]) and that equa- 
tion ([18)) shows that cooling losses should make the feedback 
ineflficient for densities nn < 31 cm"^ (T/IO'^-^)^/^ 

The wind velocities are a factor of a few higher than 
for the low-mass galaxy and, at least for AT > 10^"^ K, de- 
pend strongly on the temperature increase. Figure [3] showed 
that for these high heating temperatures, the wind blows 
channels through which it can freely escape to very large 
distances. The effect of drag by halo gas is therefore smaller 
than for the dwarf galaxy and for the low- AT versions of 
the massive galaxy, which all predict much puffier gas disks. 
Hence, in the high- AT simulations of the massive galaxy 
the outflow speed is predominantly regulated by gravity (as 
opposed to gas drag) and by the initial velocity. While the 
potential is always the same, the initial wind velocity is set 
by the heating temperature. 

The radial dependence of the outflow is shown in the 
bottom row of Fig. [9] Interestingly, the models that predict 
similar SF histories, also predict similar outflow rates for r < 
0.2rvir. However, at larger radii the outflow rates diverge, 
with higher heating temperatures giving larger mass fluxes. 
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Except for model G12-100-65, the velocities increase linearly 
with radius beyond about 0.5 rvir, suggesting that they are 
determined by travel time constraints, as was the case for 
the dwarf galaxy. 



6.2 Comparison with the kinetic feedback model 

In this section we will compare the results of simulations 
using our new thermal feedback prescription with runs em- 
ploying the kinetic feedback of DS08. As discussed in detail 
in DS08, the kinetic feedback recipe works as follows. Once 
a star particle reaches an age of 3 x 10^ yr, its neighbouring 
gas particles each have a probability of rjra^/T^^^f^ rai of re- 
ceiving a randomly oriented kick of velocity v^ . For the case 
of equal mass particles, the mass loading factor 77 equals the 
average number of particles kicked per star particle. We use 
DSOS's fiducial values of v^ = 600 kms~^ and r/ = 2, which 
correspond to 40 percent of the available energy. 

We compare the kinetic feedback runs to thermal feed- 
back simulations that use the same fraction of the available 
SN energy (i.e. fth = 0.4 as opposed to 1.0 for our fidu- 
cial model). We use a temperature increase of AT = 10^ K 
(as opposed to 10^"^ K for our fiducial model) as this is 
close to the post-shock temperature for a shock velocity of 
600 kms"\ 

To eliminate potential differences other than the feed- 
back recipe, we re-ran models ml2 and ml2nowind of DS08 
with the new code and employing the same softening lengths 
as used here. 



6.2.1 Star formation history 

Fig. [10] compares the SF histories of the thermal feedback 
runs G 10- 04 0-70 and 01 2-040-70 with the equivalent ki- 
netic feedback runs mlO and ml 2 of DS08. 

We first verify that the new criterion for identifying 
star- forming gas (see Section [5?T]) gives results that are con- 
sistent with the implementation used in DS08. In order to 
eliminate differences due to the feedback implementation, 
we compare runs without feedback. Comparison of the black 
and orange dotted curves in Fig. IIOI which show the SF his- 
tories predicted with the new and old prescriptions, respec- 
tively, shows that differences due the slight change in the 
recipe for SF are negligible. 

Comparing the orange and black solid curves, which 
show the SF histories in the runs with kinetic and ther- 
mal feedback, respectively, we see that the two methods for 
injecting the energy from SNII are generally in good agree- 
ment. 



6.2.2 Mass outflow rate and wind velocity 

Fig. [11] shows the mass outflow rate (left column) and the 
mean outflow velocity (right column) for models GIO-O4O- 
70 and mlO. The top row shows the evolution measured at 
r = 0.2rvir, while the bottom row illustrates the dependence 
on radius at time t = 500 Myr. The agreement is generally 
excellent. This implies that, at the resolution used for the 
dwarf galaxy, the average quantities of the outflow are in- 
sensitive to the form in which the energy is injected. Note 



that this is not a consequence of a fortunate choice of pa- 
rameters, as we showed in Section [6.1 .41 that the outflow is 
insensitive to the temperature increase AT. The agreement 
between simulations injecti ng kinetic and thermal energy is 
consistent with the results of lDurier fc Dalla Vecchial (|2012l l . 
Fig. [T2I shows the same plots for the 10^^/i~^ M© halo. 
While the agreement for the velocities is again very good, 
there are in this case some noticeable differences between 
the mass outflow rates. Compared with the thermal feed- 
back simulation, the kinetic model predicts a higher outflow 
rate at radii 0.2 < r/r^iv ^ 0.4 at late times. The differences 
may imply that the resolution is too low to achieve con- 
vergence between kinetic and thermal feedback prescription 
(recall that the mass resolution is two orders of magnitude 
lower than for the dwarf galaxy). Another reason, is how- 
ever, that, unlike for the dwarf galaxy, for the massive galaxy 
the outflow does depend on the wind parameters. Moreover, 
for this massive halo AT = 10^ K (or v^ = 600 kms~^) 
marks the transition between the regimes of galactic foun- 
tains (lower AT or t'w) and efficient large-scale winds (higher 
AT or Vw), as can be seen from Figs. [71 and [9l Hence, the 
results are very sensitive to small differences in the input 
parameters and we could have obtained better agreement 
by fine-tuning the value of AT (or v^^). 

6.3 Resolution tests 

We tested the numerical convergence of our implementation 
of thermal feedback by decreasing the particle numbers by 
factors of 8 and 64. We will denote the corresponding runs by 
appending 'LR008' resp. 'LR064' to the simulation names. 
Hence, the particle masses are increased by factors of 8 and 
64, while the gravitational softening lengths and, for a fixed 
density, the SPH smoothing kernels are increased by factors 
of 2 and 4, respectively. Before showing the results, it is 
useful to consider what we may expect. To do so, we have to 
check whether the simulations resolve the Jeans scales and 
whether we expect radiative cooling losses to be s ignificant. 

As discussed in lSchave fc Dalla Vecchial (|2008|), for star- 
forming gas in our fiducial simulation of the massive halo the 
ratio of the SPH kernel mass to the Jeans mass is smaller 
than 1/6 and the ratio of the SPH kernel size to the Jeans 
length is at most 1/(48) ^^^ ^ 0.28. For the low-mass halo 
the mass and length ratios are lower by factors of 100 and 
100^/"^ ~ 4.64, respectively. Note that the maximum pos- 
sible values of these ratios are independent of the density 
because star-forming particles cannot have temperatures be- 
low a power-law effective equation of state with polytropic 
index 7eff = 4/3. Thus, while even our lowest-resolution sim- 
ulation of the dwarf galaxy resolves the Jeans scales, the 
same is only true for our highest-resolution simulation of 
the massive galaxy. 

In Section H we demonstrated that, for AT = lO'^"^ K 
and our fiducial resolution for the massive galaxy, radia- 
tive losses should have little impact on the efficiency of 
the feedback if the energy is injected in gas with density 
nu < 31 cm~^ and that this critical density is inversely pro- 
portional to the squareroot of the particle mass (eq. [H]). 
Hence, for the intermediate- and low-resolution models the 
densities above which radiative losses may prevent efficient 
feedback are about 11 and 4 cm~^, respectively, and for the 
low-mass galaxy these densities are 10 times higher. Com- 
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Figure 10. Comparison of the SF histories in simulations employing thermal (GlO-040-70 and G12-040-70) and kinetic (mlO and ml2) 
feedback. All models inject 40 percent of the SNII energy. The kinetic feedback assumes an initial wind velocity of 600 kms"-*^ and the 
thermal feedback a temperature increase of 10^ K, which is close to the post-shock temperature for a shock velocity of 600 kms"-*^. The 
left and right panels show the SFR as a function of time for the 10"'^° and 10^'^h~^ Mq haloes, respectively. The very small difference 
between the no-feedback models (that are only noticeable for the massive galaxy) are due to the different treatment of star-forming gas 
(see section [5T]) . The two feedback implementations result in very similar SF histories. 
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Figure 11. Comparison of the mass outflow rate (left column) and average outflow velocity (right column) measured through a spherical 
shell at radius r = 0.2rvir as a function of time (top row) and at t = 500 Myr as a function of radius (bottom row) between model 
GlO-040-70, which employs thermal feedback, and model mlO, which uses kinetic feedback. The dotted curve in the top-left panel 
indicates the SFR of model GlO-040-70. All other curves are labelled in the legends. The agreement between the outflows predicted by 
simulations that inject the SNII energy in thermal and kinetic form is generally excellent. 



paring this to the actual densities of the star- forming gas in 
the high-resolution simulations (Figs. [Hand O, we see that 
radiative losses may not be negligible for the lower-resolution 
simulations of the high-mass galaxy, but that the feedback 
should remain efficient in the low-resolution models of the 
dwarf galaxy. 



Thus, both Jeans and radiative cooling arguments sug- 
gest that even the lowest-resolution simulation of the dwarf 
galaxy should give converged results. On the other hand, 
we do expect the lower- resolution versions of the massive 
galaxy to show some difference. Furthermore, given that the 
lowest-resolution models have fewer than 3700 gas particles 
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Figure 12. As Fig.[TT] but for the 10^'^h~^ Mq halo. The agreement between the thermal and kinetic implementations of SNII feedback 
is very good for the wind velocities, but the mass outflow rates differ substantially at late times and intermediate radii. This probably 
reflects the fact that these wind parameter values yield results that are intermediate between the galactic fountain and eflicient, large-scale 
winds regimes, making the outcome very sensitive to the exact values of the wind parameters. 
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Figure 13. Numerical convergence of the SF histories of the 10^^h~^ Mq and 10^^h~^ Mq haloes (left- and right panels, respectively). 
The mass resolution is decreased by factors of 8 (dashed curves) and 64 (dash-dotted curves). The convergence is good, although there 
is a small, systematic increase of the SFR with decreasing resolution for the massive galaxy. 



in the disk, we expect the results to become noisy due to 
the poor sampling. 

Figs. [13] shows how the SF history depends on the nu- 
merical resolution for the fiducial models GlO-100-75 (left) 
and G12-100-75 (right). Similarly, Fig. [H illustrates the res- 
olution dependence of the predicted evolution of the mass 
outflow rate (left column) and wind velocity (right column) 
for the dwarf (top row) and massive (bottom row) galaxies. 
Clearly, our expectations are borne out. 

Whereas the dwarf galaxy runs are very well converged. 



for the massive galaxy the predicted SFR increases with 
resolution, although the effect is small. The first 250 Myr the 
mass outflow rate is somewhat higher in the lower resolution 
simulations of the massive galaxy, but the situation reverses 
at later times. Except for the first 10 Myr, the differences 
are, however, small. The wind velocity is more sensitive to 
the resolution and is generally higher in the lower-resolution 
models. The predictions for the outflow properties become 
noisy for the lowest resolution simulations. 

It is interesting that the convergence with numerical 
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Figure 14. Resolution dependence of tlie mass outflow rate (left column) and mean outflow velocity (right column) measured through a 
spherical shell at radius r = 0.2rvir as a function of time for the 10-'^'-' (top row) and the 10^^h~^ Mq (bottom row) haloes. Note that the 
particle mass in G10-100-75-LR064 is still lower than that in G12-100-75. While the predictions for the low-mass galaxy are converged, 
decreasing the resolution yields mostly higher outflow velocities for the high-mass galaxy, but has little effect on the mass outflow rates. 



resolution is better than found by DS08 for the kinetic feed- 
back versions of the same simulations (c.f. Figs. 10 and 11 
of DS08). For example, with kinetic feedback the outflow 
rate for the massive galaxy increased by a factor of about 
3 as the mass resolution was decreased by a factor of 8 and 
the SFR of the dwarf galaxy increased substantially when 
the particle mass was decreased by a factor of 64. It should 
be kept in mind, however, that this could in part be due 
to the fact that the fiducial wind velocity used by DS08 
(600 kms""*^) corresponds to post-shock temperature jumps 
that are a factor of a few lower than our fiducial value of 
AT = 10^-^ K. 



7 DISCUSSION 

Models of galaxy formation and evolution require feedback 
from star formation to reproduce the observed properties of 
galaxies. Cosmological simulations do not have sufficient res- 
olution to resolve individual SN explosions and must there- 
fore resort to sub-grid recipes. The simplest method, inject- 
ing the SN energy released by a star particle (i.e. a simple 
stellar population; SSP) during each time step into its sur- 
roundings, does not lead to efficient feedback because the 
injected thermal energy is quickly radiated away. Success- 
ful recipes generally either inject the energy in kinetic form, 
turn off radiative cooling temporarily, or inject the energy 
in a hot sub-grid phase that is decoupled from the colder 
phases by hand. 



We demonstrated that the catastrophic radiative losses 
suffered by simple thermal feedback recipes are due to a 
mismatch between the gas mass in which the energy is 
injected and the mass of the SSP that produced the en- 
ergy. In the real Universe, one SNII is produced for ev- 
ery r^ 10^ M© of stars, and the energy released in the 
explosion is initially carried by ^ 10^ M© of eject a. Be- 
cause the ratio between the mass of the eject a and the 
mass of the stellar population that released the energy is 
small (^ 1), the SN energy per unit mass of ejecta is 
high. Hence, the ejecta move at very high velocity (^ 
10^ kms~^) and the post-shock temperatures are suffi- 
ciently large for the radiative cooling time to be long. In- 
deed, observations do provide evi dence for hot gas associated 
with f ast galactic outflows (e. g. iHeckman, Armus, k, Milevl 
1987; Strickland et al. 20 0Q|; see a lso the review by 
Veilleux, Cecil, & Bland-Ha wthorn| [20Q5) . 

In simulations, however, the ratio between the gas mass 
receiving the SN energy and the stellar mass that produced 
it is much larger. Even if all the SNII energy of an SSP is 
injected at once, the ratio between the mass of the 'ejecta' 
and that of the SSP will typically be large. For example, 
if, for the case of SPH, the energy is shared by all A^sph 
neighbours, then the ratio will be A^sph ^ 10 and even 
larger if multiple star particles are spawned per gas parti- 
cle or if the feedback energy is released over multiple time 
steps. Consequently, the temperature of the heated particles 
will be relatively low and most of the thermal energy will 
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be radiated away before it is able to do PdV work on its 
surroundings. It is important to note that the ratio between 
the gas mass receiving the energy and the mass of the SSP 
that produced it is independent of the resolution. 

The realisation that the cause of the inefficient thermal 
feedback is a mismatch between the simulated and observed 
ratios of heated mass to stellar mass, rather than a straight- 
forward lack of numerical resolution, immediately indicates 
the solution: we need to decrease this ratio in the simula- 
tions. By decreasing the mass ratio, we increase the temper- 
ature of the heated gas and hence its cooling time. If the 
cooling time is long compared with the sound crossing time 
scale across a resolution element, the heated gas will begin 
to expand adiabatically and the injected thermal energy will 
be efficiently converted into kinetic energy. We showed that 
in the temperature regime for which Brehmsstrahlung dom- 
inates the radiative cooling rate ( ^ 10^ K for solar abun- 
dances) , the ratio of the cooling time and the sound crossing 
time remains constant if the gas expands adiabatically, so 
that adiabatic cooling does not invalidate the argument. 

We showed analytically (see eq. [18]) that, for the case 
of SPH, the maximum density for which radiative losses 
are small is na - 26 cm-^{T/10^-^ K)^^^{m/10^ M0)-^/^ 
where T is the temperature of the heated gas after receiving 
the feedback energy and m is the mass of a gas resolution 
element (and we assumed that AT is sufficiently high for 
Brehmsstrahlung to dominate the radiative cooling). The 
maximum density is nearly the same for AMR simulations 
with cell sizes that are at least a factor of 4 smaller than the 
local Jeans length (evaluated in gas with this density and 
a temperature of 10^ K, see eq. [20]). Hence, by specifying 
the desired temperature jump AT, we can guarantee that 
the feedback is efficient up to some gas density that we can 
estimate analytically. 

We implemented this idea in the SPH code gadget in 
the form of stochastic thermal feedback. A fixed AT then 
translates into a fixed probability of receiving energy, which 
we evaluate for each SPH neighbour of a star particle that 
has just crossed the critical age tsN = 3x 10^ yr, correspond- 
ing to the maximum lifetime of stars that end their lives as 
SNII. For a Chabrier IMF and assuming equal mass parti- 
cles, the expectation value of the number of heated particles 
is 1.34/th(AT/10'^-^ K)-\ where fth is the fraction of the 
SNII energy that is injected (see eq. [8]). 

Note that the parameter AT plays a similar role as the 
initial wind velocity Vw for the case of kinetic feedback. For 
a fixed I'w, the mass loading factor used in kinetic feedback 
implementations sets fth- The fraction of the injected SN 
energy can also be used as a second free parameter for the 
case of thermal feedback. For a fixed temperature jump, the 
"initial mass loading" is then the ratio of the heated mass 
per unit stellar mass formed and this ratio is proportional 
to the parameter fth- 

The combination of stochastic feedback and a fixed in- 
crease in the energy per unit mass of the gas receiving the 
feedback energy also works for kinetic feedback: we can spec- 
ify Vw and give each neighbouring resolution element of a 
star particle a probability of being kicked in the wind that 
is proportional to the mass loading factor or, equivalent ly, 
to fth - This is in fact exactly the implementation of kinetic 
feedback that we used in DS08. 

The equivalence of the roles of AT and t'w suggests 



that the efficiency of the kinetic feedback does not depend 
directly on the ratio of the initial wind velocity and the 
escape velocity, contrary to what is often assumed. Indeed, 
in DS08 we showed that if t'w is too low, the wind stalls 
in the ISM, before it has even begun to climb out of the 
gravitational potential. As we have shown, the efficiency of 
the wind depends on the radiative losses and hence, for a 
fixed value of t'w, on the numerical resolution and the gas 
density. Because the typical gas densities increase with the 
gas pressure and thus with the depth of the potential well, 
the escape velocity does matter indirectly (also because drag 
forces increase with the pressure). If the wind manages to 
blow out of the ISM, then the ultimate efficiency of the 
feedback does depend on the ratio of the velocity of the 
wind leaving the ISM and the escape velocity, because the 
ejected gas will rain back onto the galaxy if it cannot escape 
the galaxy's potential well. 

We presented analytic derivations of the resolution cri- 
teria, both for SPH and AMR simulations. We tested our 
recipe for thermal feedback on SPH simulations of isolated 
disc galaxies in dark matter haloes of total mass 10h~^ M© 
and 10^^ h~^ M© using the same set-up as we used to study 
kinetic SN feedback in DS08. We explored the effect of the 
feedback on the gas distribution, the star formation history, 
the mass outflow rate, and the wind velocity. The results 
were in accord with our analytic predictions. 

For sufficiently high AT and for sufficiently high resolu- 
tion, the thermal feedback strongly reduces the star forma- 
tion rate and results in a strong, large-scale, bi-polar out- 
flow. Reassuringly, the results converge with both AT and 
resolution and the converged results also agree well with sim- 
ulations employing kinetic feedbacHJ (with sufficiently high 

However, if AT and/or the resolution are too low, then 
the results become sensitive to both. For a fixed resolution, 
higher values of AT will then result in more efficient feed- 
back. Hence, in this regime the ability to choose AT or, for 
the case of kinetic feedback t'w , implies a considerable free- 
dom. This freedom associated with the implementation of 
feedback from star formation is currently the limiting fac- 
tor for the predictive power of cosmological simulations (e.g. 
ISchave et aIll2010l : [ScannaDieco et al.ll2Qllh . 

Given that a higher AT yields smaller radiative losses, 
one may ask why we do not use ultra- high values. Indeed, 
even in low-resolution simulations the feedback could be 
made efficient locally by increasing AT (or v^^ for the case of 
kinetic feedback). There are, however, several reasons why it 
is undesirable to increase this parameter to values > 10^ K. 
First, even for fth = 1 such high heating temperatures im- 
ply that, on average, each star particle will heat less than 
one neighbouring gas particle. Such a situation breaks the 
locality of the feedback and may lead to sampling problems. 
That this can have grave consequences is easy to see by con- 
sidering the limiting case in which the number of heated gas 
particles per star particle (i.e. the mean initial mass loading) 
is ^ 1. Most heavy elements released by massive stars will 



^ The thermal feedback models only agree with kinetic feedback 
simulations if the kicked wind particles are not temporarily decou- 
pled from the hydrodynamics, see DS08 for a critical discussion 
of this common practice. 
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then no longer be injected in a wind. Many generations of 
star particles can form in a given gas cloud before a single 
feedback event takes place. Conversely, if, despite the low 
probability, a star particle forming in a region with a low 
star formation density does heat a neighbour, the energy in- 
jected may be sufficiently large to do catastrophic damage. 
Clearly, we should avoid the regime in which the expecta- 
tion value for the number of heated gas elements per star 
particle formed is much less than one, at least for galaxies 
resolved with relatively small numbers of particles. 

At present, large- volume cosmological simulations typi- 
cally have particle masses ^ 10^ M0. At this resolution even 
heating one gas particle per star particle (which corresponds 
to AT - 10'^-^ K for a Chabrier IMF and fth = 1) results in 
strong radiative losses for densities nn > 10 cm~^, which are 
routinely reached in such simulations. Hence, the predictions 
are still sensitive to the values of the feedback prescription 
and are thus uncertain. This undesirable limitation can be 
turned into an advantage if one takes an approach similar in 
spirit to semi-analytic models: by varying AT (or v^r) with 
halo mass or with the local physical conditions, the feedback 
can be tuned to reproduce the desired galaxy formation ef- 
ficiency. However, the arguments given above demonstrate 
that the results may change with increasing resolutiorlj. If 
the resolution and the value of AT are sufficiently high for 
cooling losses to be small, then the feedback can still be 
tuned by varying fth with the local conditions. 

We have demonstrated that, contrary to common wis- 
dom, thermal feedback can be efficient without turning off 
radiative cooling. For sensible parameter choices overcooling 
can already be avoided for densities typical of the warm ISM 
(i.e. nu ^ 1 cm~^) at the resolution achievable for large- 
scale cosmological simulations (m ^ 10^ M0) and for sim- 
ulations of individual (low- to intermediate-mass) galaxies 
we can already afford the resolution (m ^ 10^ M©) required 
for the feedback to remain efficient up to densities typical of 
molecular clouds (nu '^ 10^ cm""^). We have also shown that 
with sufficient resolution, the results become insensitive to 
the problematic parameter of the feedback implementation 
(i.e. AT for thermal feedback and v^ for kinetic feedback) 
and the form in which the energy is injected, thus removing 
some of the most important uncertainties in the ingredients 
of hydrodynamical simulations of galaxy formation. 



ACKNOWLEDGEMENTS 

We are very grateful to Volker Springel for allowing us to use 
GADGET and his initial conditions code for the simulations 
presented here. We thank Rob Crain and Jarrett Johnson 
for a careful reading of the manuscript, and the anonymous 
referee for a helpful report. The simulations presented here 
were run on the Cosmology Machine at the Institute for 
Computational Cosmology in Durham as part of the Virgo 
Consortium research programme and on the TMoX cluster 



^ Turning off the hydrodynamical forces in the high-density 
regime co uld make the result s insensitive to resolution for kinetic 
feedback (jSpringel fc Hernauist 2003) . but the predictions will in 
that case disagree with converged, self-consistent high-resolution 
simulations. 



at the Rechenzentrum Garching of the Max Planck Soci- 
ety. This work was supported by Marie Curie Reintegration 
Grant FP7-RG-256573 and by the Marie Curie Initial Train- 
ing Network CosmoComp (PITN-GA-2009-238356). 



REFERENCES 

Booth C. M., Schaye J., 2009, MNRAS, 398, 53 

Brook C. B., Kawata D., Gibson B. K., Flynn C, 2004, 

MNRAS, 349, 52 
Chiosi C, Bertelh G., Bressan A., 1992, ARA&A, 30, 235 
Creasey P., Theuns T., Bower R. G., Lacey C. G., 2011, 

MNRAS, 415, 3706 
Dalla Vecchia C, Schaye J., 2008, MNRAS, 387, 1431 

(DS08) 
Dubois Y., Teyssier R., 2008, A&A, 477, 79 
Durier F., Dalla Vecchia C, 2012, MNRAS, 419, 465 
Ferland, G. J. 2000, Revista Mexicana de Astronomia y 

Astrofisica Conference Series, 9, 153 
Gerritsen J. P. E., 1997, PhDT 
Haardt, F., & Madau, P. 2001, in the proceedings of 

XXXVI Rencontres de Moriond, |astro-p h/0106018' 
Heckman T. M., Armus L., Miley G. K., 1987, AJ, 93, 276 
Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 

421, 3522 
Kawata D., 2001, ApJ, 558, 598 
Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 

19 
Kay S. T., Pearce F. R., Frenk C. S., Jenkins A., 2002, 

MNRAS, 330, 113 
Kay S. T., Thomas P. A., Theuns T., 2003, MNRAS, 343, 

608 
Kennicutt R. C, Jr., 1998, ApJ, 498, 541 
Marri S., White S. D. M., 2003, MNRAS, 345, 561 
Mihos J. C, Hernquist L., 1994, ApJ, 437, 611 
Monaghan J. J., 1997, JCoPh, 136, 298 
Mori M., Yoshii Y., Tsujimoto T., Nomoto K., 1997, ApJ, 

478, L21 
Murante G., Monaco P., Giovalli M., Borgani S., Diaferio 

A., 2010, MNRAS, 405, 1491 
Navarro J. F., White S. D. M., 1993, MNRAS, 265, 271 
Oppenheimer B. D., Dave R., 2006, MNRAS, 373, 1265 
Osterbrock D. E., 1989, agna.book, 
Piontek F., Steinmetz M., 2011, MNRAS, 410, 2625 
Saitoh T. R., Makino J., 2009, ApJ, 697, L99 
Scannapieco C, Tissera P. B., White S. D. M., Springel V., 

2006, MNRAS, 371, 1125 
Scannapieco C, Wadepuhl M., Parry O. H., et al., 2011, 

larXiv:1112.0315 l 
Schaye J., 2004, ApJ, 609, 667 

Schaye J., Dalla Vecchia C, 2008, MNRAS, 383, 1210 
Schaye J., Dalla Vecchia C, Booth C. M., et al., 2010, 

MNRAS, 402, 1536 
Sommer-Larsen J., Gotz M., Portinari L., 2003, ApJ, 596, 

47 
Springel V., Hernquist L., 2003, MNRAS, 339, 289 (SH03) 
Springel V., 2005, MNRAS, 364, 1105 
Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 

361, 776 
Stinson G., Seth A., Katz N., Wadsley J., Governato F., 

Quinn T., 2006, MNRAS, 373, 1074 



© 0000 RAS, MNRAS 000, 000-000 



Galactic outflows with thermal SN feedback 21 



Strickland D. K., Heckman T. M., Weaver K. A., Dahlem 

M., 2000, AJ, 120, 2965 
Thacker R. J., Couchman H. M. P., 2000, ApJ, 545, 728 
Theuns T., Leonard A., Efstathiou G., Pearce F. R., 

Thomas P. A., 1998, MNRAS, 301, 478 
Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 

43, 769 
Wiersma R. P. C., Schaye J., Smith B. D., 2009a, MNRAS, 

393, 99 
Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., 

Tornatore L., 2009b, MNRAS, 399, 574 



© 0000 RAS, MNRAS 000, 000-000 



